Prepare various bits of data

This notebook prepares the - locational data for patents - person/firm level identifiers

Load data

condensed address data…

We are working out the exact location for a sample of firms. This exploits cheap/free usage of APIs. Specifically, the mapbox API allows 100k free calls per month. We use the exactly matched address strings to assign an approximate address to the remaining persons in the patent data. Specifically, we assign a random location of the found addresses in the same country to all addresses with country information. Note that this could be greatly refined in future work; e.g. predict the likely location use features such as sector etc. The main purpose of the exercise is now to get map points that are reasonable within specific countries.

get address data from bigquery unless already stored locally….

personfile=paste0(localbig,"\\data\\person_address.parquet")
if (!file.exists(personfile)) {
  library(bigrquery)
  library(DBI)
  
  # Set your project ID
  project_id <- "patbis"
  
  # Option 1: Direct table download
  # Specify dataset and table name
  dataset <- "inglobe"
  table <- "person_address_extract"
  
  # Download the entire table
  df <- bq_table_download(
    bq_table(project_id, dataset, table)
  )
  
  
  library(arrow)
  write_parquet(df,personfile)
}

Get person data

a little check…. are high value patents MNE ones?

#bigX=countrymap %>% inner_join(hicinno) %>% inner_join(istraxglobal)


bigX=countrymap %>% inner_join(hicinno) %>% 
   inner_join(istraxglobal) %>%
   inner_join(istraxemdeexcn) %>% 
   inner_join(istraxnational) %>% 
   rename(iso2=ctry_code) %>% 
   inner_join(country_income)%>% mutate(CNid=iso2=="CN") %>% left_join(greenid) %>% 
   mutate(anygreen=!is.na(technology))

Joining with `by = join_by(docdb_family_id)`
Joining with `by = join_by(docdb_family_id, ctry_code)`
Joining with `by = join_by(docdb_family_id, ctry_code)`
Joining with `by = join_by(docdb_family_id, ctry_code)`
Joining with `by = join_by(iso2)`
Joining with `by = join_by(docdb_family_id)`

pie charts

# Create summary
collab_summary <- bigX %>%  filter(HIC==0&CNid==0 & !(hicinno == 1 & multi == 0)) %>% 
  count(hicinno, multi) %>%
  mutate(category = case_when(
    hicinno == 1 & multi == 1 ~ "HIC & International",
    hicinno == 1 & multi == 0 ~ "HIC Only",
    hicinno == 0 & multi == 1 ~ "International Only",
    hicinno == 0 & multi == 0 ~ "No Collaboration"
  ))

# Create interactive pie chart
plot_ly(collab_summary, 
        labels = ~category, 
        values = ~n, 
        type = 'pie',
        textinfo = 'label+percent+value',
        hoverinfo = 'text',
        text = ~paste('Count:', n)) %>%
  layout(title = 'Distribution of Collaboration Types')

donut

simple regressions

model=feols(istrax_EMDEexCN*100 ~ hicinno+multi,
      bigX %>% filter(HIC==0&CNid==0))

print(model)

OLS estimation, Dep. Var.: istrax_EMDEexCN * 100
Observations: 323,210
Standard-errors: IID 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 46.4   Adj. R2: 0.006315

plot_collab_returns <- function(data = bigX, 
                                var, 
                                filter_HIC = TRUE,
                                filter_CNid = TRUE,
                                exclude_HIC_only = TRUE,
                                title = "Return by Collaboration Type",
                                xlab = "Spillover return per $100 invested",
                                multiply_by = 100,
                                palette = "Set1",
                                text_size = 3.5,
                                label_size = 4) {
  
  var_enquo <- enquo(var)
  
  # Apply filters conditionally
  filtered_data <- data
  if (filter_HIC) filtered_data <- filtered_data %>% filter(HIC == 0)
  if (filter_CNid) filtered_data <- filtered_data %>% filter(CNid == 0)
  if (exclude_HIC_only) filtered_data <- filtered_data %>% filter(!(hicinno == 1 & multi == 0))
  
  # Prepare summary data
  summary_data <- filtered_data %>%
    mutate(collab_type = case_when(
      hicinno == 1 & multi == 1 ~ "HIC & International",
      hicinno == 1 & multi == 0 ~ "HIC Only",
      hicinno == 0 & multi == 1 ~ "International Only",
      hicinno == 0 & multi == 0 ~ "No Collaboration"
    )) %>%
    group_by(collab_type) %>%
    summarise(
      avg_value = mean(!!var_enquo * multiply_by, na.rm = TRUE),
      n = n(),
      nn = log(n),
      .groups = "drop"
    ) %>%
    mutate(
      share = n / sum(n),
      log_share = nn / sum(nn),
      log_share_scaled = (log_share - min(log_share)) * 50,
      y_pos = row_number(),
      ymin = y_pos - log_share/2,
      ymax = y_pos + log_share/2
    ) %>%
    arrange(desc(avg_value))
  
  # Create plot
  ggplot(summary_data) +
    geom_rect(aes(xmin = 0, xmax = avg_value, ymin = ymin, ymax = ymax, fill = collab_type),
              alpha = 0.8, color = "white", linewidth = 1) +
    geom_text(aes(x = avg_value + max(avg_value)*0.02, y = y_pos, 
                  label = paste0(round(avg_value, 2))),
              hjust = 0, size = label_size, fontface = "bold", color = "black") +
    geom_text(aes(x = 0, y = y_pos, 
                  label = paste0(collab_type, "\n(n=", scales::comma(n), ", ", 
                                round(share*100, 1), "%)")),
              hjust = 1.05, size = text_size, fontface = "bold", color = "black") +
    scale_fill_brewer(palette = palette) +
    labs(title = title,
         subtitle = "Bar width proportional to log(share) of observations",
         x = xlab,
         y = NULL) +
    theme_minimal() +
    theme(
      legend.position = "none",
      panel.grid.major.y = element_blank(),
      plot.title = element_text(face = "bold", size = 14),
      axis.text.y = element_blank(),
      plot.margin = margin(10, 10, 10, 100)
    ) +
    coord_cartesian(clip = "off")
}

# Usage:
plot_collab_returns(bigX, istrax_global)

plot_collab_returns(bigX, istrax_EMDEexCN)

NA
NA
NA

find locations…

we are using tidygeocoder & mapbox…

library(arrow)
df=read_parquet(paste0(localbig,"\\data\\person_address.parquet"))

library(tidygeocoder)

# create a dataframe with addresses

## let's check if we have stored such a dataframe previously
mboxfile=paste0(localbig,"\\data\\lat_longs_mapbox.parquet")
if (!file.exists(mboxfile)) {
  #df=df %>% mutate(rr=runif(1:n()))  %>% arrange(rr) %>% select(-rr)
  some_addresses=df %>% head(10) %>% mutate(addressstr=paste(person_address,person_ctry_code,sep=","))
  # geocode the addresses
  lat_longs_mapbox <- some_addresses %>%
    geocode(addressstr, method = 'mapbox', 
            #flatten=T,
            #full_results=T,
            #id = id,
            
            #type= type,
            lat = latitude , long = longitude)
  #> Passing 3 addresses to the Nominatim single address geocoder
  #> Query completed in: 3.2 seconds
} else{
  
  lat_longs_mapbox=read_parquet(mboxfile)
  
}



recursive_mapbox=function(){
  df_filtered <- df %>%
    anti_join(lat_longs_mapbox, by = "person_id")
  
  df_filtered=df_filtered %>% mutate(rr=runif(1:n()))  %>% arrange(rr) %>% select(-rr)
  
  some_addresses=df_filtered %>% head(1500) %>%  mutate(addressstr=paste(person_address,person_ctry_code,sep=","))
  print(nrow(df_filtered))
  lat_longs_mapbox2 <- some_addresses %>%
    geocode(addressstr, method = 'mapbox', lat = latitude , long = longitude,full_results=F)
  
  
  lat_longs_mapbox_new=lat_longs_mapbox %>% bind_rows(lat_longs_mapbox2)
  #lat_longs_mapbox=lat_longs_mapbox %>% bind_rows(lat_longs_mapbox2)
  return(lat_longs_mapbox_new)
}

#lat_longs_mapbox=lat_longs_mapbox_new  %>% distinct()
#lat_longs_mapbox=lat_longs_mapbox  %>% distinct() %>% select(-rr)



lat_longs_mapbox=recursive_mapbox()
write_parquet(lat_longs_mapbox,paste0(localbig,"\\data\\lat_longs_mapbox.parquet"))

reading back in from disk

lat_longs_mapbox=read_parquet(paste0(localbig,"\\data\\lat_longs_mapbox.parquet"))

Cleaning up

We examine if the points suggested fall within the country boundaries, which doesn’t always seem the case…

test=lat_longs_mapbox %>% filter(nchar(person_address)<=2)

testdf=lat_longs_mapbox %>% filter(!is.na(person_ctry_code) & !is.na(longitude))


library(sf)
library(rnaturalearth)

# Get world boundaries
world <- ne_countries(scale = "medium", returnclass = "sf")








world <- world %>%
  mutate(iso_a2_fixed = case_when(
    name == "Taiwan" ~ "TW",
    name == "Kosovo" ~ "XK",
    iso_a2 == "-99" & !is.na(iso_a3) ~ iso_a3,  # Fallback to iso_a3
    TRUE ~ iso_a2
  ))










# Convert dataframe to spatial points
df_sf <- st_as_sf(testdf, 
                  coords = c("longitude", "latitude"), 
                  crs = 4326)

prepare progress bar

library(progress)

n <- nrow(testdf)

pb <- progress_bar$new( format = “[:bar] :percent eta: :eta”, total = n, clear = FALSE, width = 60 )

Initialize ok column

testdf$ok <- FALSE

Check each point

for(i in 1:nrow(testdf)) { country_poly <- world[world\(iso_a2 == testdf\)person_ctry_code[i], ] if(nrow(country_poly) > 0) { testdf\(ok[i] <- st_intersects(df_sf[i,], country_poly, sparse = FALSE)[1,1] } pb\)tick() }

# Detect number of cores
library(parallel)
n_cores <- detectCores() - 2  # Leave one core free

# Create cluster
cl <- makeCluster(n_cores)

# Export necessary objects to cluster
clusterExport(cl, c("df_sf", "world", "testdf"), envir = environment())
clusterEvalQ(cl, library(sf))

# Parallel processing
testdf$ok <- parSapply(cl, 1:nrow(testdf), function(i) {
  country_poly <- world[world$iso_a2 == testdf$person_ctry_code[i], ]
  if(nrow(country_poly) > 0) {
    return(st_intersects(df_sf[i,], country_poly, sparse = FALSE)[1,1])
  } else {
    return(FALSE)
  }
})


# Stop cluster
stopCluster(cl)


test=testdf %>% filter(person_ctry_code=="TW")


# correct some issues
testdf=testdf %>% mutate(ok=ifelse(person_ctry_code=="TW" ,T,ok), # correct some Chines corruption of the shapefile
                         ok=ifelse(nchar(gsub("\\s+", "", person_ctry_code))<2 ,T,ok) )

Infer iso2

library(sf)
library(rnaturalearth)
library(dplyr)

# Get world boundaries
world <- ne_countries(scale = "medium", returnclass = "sf")

# Initialize inferrediso2 column
lat_longs_mapbox <- lat_longs_mapbox %>%
  mutate(inferrediso2 = NA_character_)

# Filter to only rows with valid coordinates
valid_coords <- !is.na(lat_longs_mapbox$latitude) & 
                !is.na(lat_longs_mapbox$longitude)

if(sum(valid_coords) > 0) {
  # Convert to spatial points (only valid coordinates)
  points_sf <- st_as_sf(lat_longs_mapbox[valid_coords, ], 
                        coords = c("longitude", "latitude"), 
                        crs = 4326)
  
  # Spatial join to find which country polygon contains each point
  points_with_country <- st_join(points_sf, 
                                 world %>% select(iso_a2), 
                                 join = st_intersects,
                                 left = TRUE)
  
  # Extract inferrediso2
  lat_longs_mapbox$inferrediso2[valid_coords] <- points_with_country$iso_a2
  
  # Special handling for Taiwan
  # Define Taiwan's bounding box with some buffer
  taiwan_bounds <- list(
    lat_min = 21.8,
    lat_max = 25.4,
    lon_min = 119.3,
    lon_max = 122.1
  )
  
  # Identify points in Taiwan region
  in_taiwan_region <- valid_coords & 
                      lat_longs_mapbox$latitude >= taiwan_bounds$lat_min & 
                      lat_longs_mapbox$latitude <= taiwan_bounds$lat_max & 
                      lat_longs_mapbox$longitude >= taiwan_bounds$lon_min & 
                      lat_longs_mapbox$longitude <= taiwan_bounds$lon_max
  
  # Override with TW for points in Taiwan region (regardless of what spatial join said)
  lat_longs_mapbox$inferrediso2[in_taiwan_region] <- "TW"
}

clean up some more…

lat_longs_mapbox$inferrediso2[lat_longs_mapbox$inferrediso2=="-99"|lat_longs_mapbox$inferrediso2==""] <- NA

lat_longs_mapbox= lat_longs_mapbox %>% mutate(person_ctry_code=gsub("\\s+", "", person_ctry_code)) %>%
                                       mutate(person_ctry_code=ifelse(person_ctry_code=="" & !is.na(inferrediso2),inferrediso2,person_ctry_code)) 


test=lat_longs_mapbox %>% filter(person_ctry_code=="TW")



lat_longs_mapbox=lat_longs_mapbox %>% filter(!is.na(latitude) ) %>% 
                                      filter(!is.na(person_ctry_code)) %>% 
                                      filter(person_ctry_code!="") %>% 
                                      filter(!(person_ctry_code!=inferrediso2 & !is.na(inferrediso2) ))



lat_longs_bycountry=lat_longs_mapbox %>% group_by(person_ctry_code,longitude,latitude) %>% summarise(innos=n())

Let’s have a look:

library(leaflet)

# With clustering (much better for 100k points!)
leaflet(lat_longs_mapbox) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~longitude,
    lat = ~latitude,
    radius = 3,
    fillOpacity = 0.5,
    stroke = FALSE,
    clusterOptions = markerClusterOptions()
  )

deal with capital cities …. gonna be used as default…

# Create a comprehensive dataframe of capital cities with coordinates
capital_coords <- data.frame(
  iso2 = c(
    "AF", "AL", "DZ", "AD", "AO", "AG", "AR", "AM", "AU", "AT",
    "AZ", "BS", "BH", "BD", "BB", "BY", "BE", "BZ", "BJ", "BT",
    "BO", "BA", "BW", "BR", "BN", "BG", "BF", "BI", "KH", "CM",
    "CA", "CV", "CF", "TD", "CL", "CN", "CO", "KM", "CG", "CD",
    "CR", "CI", "HR", "CU", "CY", "CZ", "DK", "DJ", "DM", "DO",
    "EC", "EG", "SV", "GQ", "ER", "EE", "ET", "FJ", "FI", "FR",
    "GA", "GM", "GE", "DE", "GH", "GR", "GD", "GT", "GN", "GW",
    "GY", "HT", "HN", "HU", "IS", "IN", "ID", "IR", "IQ", "IE",
    "IL", "IT", "JM", "JP", "JO", "KZ", "KE", "KI", "KP", "KR",
    "KW", "KG", "LA", "LV", "LB", "LS", "LR", "LY", "LI", "LT",
    "LU", "MK", "MG", "MW", "MY", "MV", "ML", "MT", "MH", "MR",
    "MU", "MX", "FM", "MD", "MC", "MN", "ME", "MA", "MZ", "MM",
    "NA", "NR", "NP", "NL", "NZ", "NI", "NE", "NG", "NO", "OM",
    "PK", "PW", "PS", "PA", "PG", "PY", "PE", "PH", "PL", "PT",
    "QA", "RO", "RU", "RW", "KN", "LC", "VC", "WS", "SM", "ST",
    "SA", "SN", "RS", "SC", "SL", "SG", "SK", "SI", "SB", "SO",
    "ZA", "SS", "ES", "LK", "SD", "SR", "SZ", "SE", "CH", "SY",
    "TW", "TJ", "TZ", "TH", "TL", "TG", "TO", "TT", "TN", "TR",
    "TM", "TV", "UG", "UA", "AE", "GB", "US", "UY", "UZ", "VU",
    "VA", "VE", "VN", "YE", "ZM", "ZW"
  ),
  capital = c(
    "Kabul", "Tirana", "Algiers", "Andorra la Vella", "Luanda", "Saint John's", "Buenos Aires", "Yerevan", "Canberra", "Vienna",
    "Baku", "Nassau", "Manama", "Dhaka", "Bridgetown", "Minsk", "Brussels", "Belmopan", "Porto-Novo", "Thimphu",
    "Sucre", "Sarajevo", "Gaborone", "Brasília", "Bandar Seri Begawan", "Sofia", "Ouagadougou", "Gitega", "Phnom Penh", "Yaoundé",
    "Ottawa", "Praia", "Bangui", "N'Djamena", "Santiago", "Beijing", "Bogotá", "Moroni", "Brazzaville", "Kinshasa",
    "San José", "Yamoussoukro", "Zagreb", "Havana", "Nicosia", "Prague", "Copenhagen", "Djibouti", "Roseau", "Santo Domingo",
    "Quito", "Cairo", "San Salvador", "Malabo", "Asmara", "Tallinn", "Addis Ababa", "Suva", "Helsinki", "Paris",
    "Libreville", "Banjul", "Tbilisi", "Berlin", "Accra", "Athens", "St. George's", "Guatemala City", "Conakry", "Bissau",
    "Georgetown", "Port-au-Prince", "Tegucigalpa", "Budapest", "Reykjavik", "New Delhi", "Jakarta", "Tehran", "Baghdad", "Dublin",
    "Jerusalem", "Rome", "Kingston", "Tokyo", "Amman", "Nur-Sultan", "Nairobi", "Tarawa", "Pyongyang", "Seoul",
    "Kuwait City", "Bishkek", "Vientiane", "Riga", "Beirut", "Maseru", "Monrovia", "Tripoli", "Vaduz", "Vilnius",
    "Luxembourg", "Skopje", "Antananarivo", "Lilongwe", "Kuala Lumpur", "Malé", "Bamako", "Valletta", "Majuro", "Nouakchott",
    "Port Louis", "Mexico City", "Palikir", "Chișinău", "Monaco", "Ulaanbaatar", "Podgorica", "Rabat", "Maputo", "Naypyidaw",
    "Windhoek", "Yaren", "Kathmandu", "Amsterdam", "Wellington", "Managua", "Niamey", "Abuja", "Oslo", "Muscat",
    "Islamabad", "Ngerulmud", "Ramallah", "Panama City", "Port Moresby", "Asunción", "Lima", "Manila", "Warsaw", "Lisbon",
    "Doha", "Bucharest", "Moscow", "Kigali", "Basseterre", "Castries", "Kingstown", "Apia", "San Marino", "São Tomé",
    "Riyadh", "Dakar", "Belgrade", "Victoria", "Freetown", "Singapore", "Bratislava", "Ljubljana", "Honiara", "Mogadishu",
    "Pretoria", "Juba", "Madrid", "Colombo", "Khartoum", "Paramaribo", "Mbabane", "Stockholm", "Bern", "Damascus",
    "Taipei", "Dushanbe", "Dodoma", "Bangkok", "Dili", "Lomé", "Nuku'alofa", "Port of Spain", "Tunis", "Ankara",
    "Ashgabat", "Funafuti", "Kampala", "Kyiv", "Abu Dhabi", "London", "Washington, D.C.", "Montevideo", "Tashkent", "Port Vila",
    "Vatican City", "Caracas", "Hanoi", "Sana'a", "Lusaka", "Harare"
  ),
  latitude = c(
    34.5553, 41.3275, 36.7538, 42.5063, -8.8383, 17.1175, -34.6037, 40.1792, -35.2809, 48.2082,
    40.4093, 25.0343, 26.0667, 23.8103, 13.0969, 53.9045, 50.8503, 17.2510, 6.4969, 27.4728,
    -19.0196, 43.8564, -24.6282, -15.8267, 4.8895, 42.6977, 12.3714, -3.3731, 11.5564, 3.8480,
    45.4215, 14.9177, 4.3947, 12.1348, -33.4489, 39.9042, 4.7110, -11.7172, -4.2634, -4.4419,
    9.9281, 6.8270, 45.8150, 23.1136, 35.1856, 50.0755, 55.6761, 11.8251, 15.3017, 18.4861,
    -0.1807, 30.0444, 13.6929, 3.7504, 15.3229, 59.4370, 9.0320, -18.1248, 60.1699, 48.8566,
    0.4162, 13.4549, 41.7151, 52.5200, 5.6037, 37.9838, 12.0561, 14.6349, 9.6412, 11.8637,
    6.8013, 18.5944, 14.0723, 47.4979, 64.1466, 28.6139, -6.2088, 35.6892, 33.3152, 53.3498,
    31.7683, 41.9028, 17.9714, 35.6762, 31.9454, 51.1694, -1.2921, 1.3382, 39.0392, 37.5665,
    29.3759, 42.8746, 17.9757, 56.9496, 33.8938, -29.3167, 6.3156, 32.8872, 47.1410, 54.6872,
    49.6116, 41.9973, -18.8792, -13.9626, 3.1390, 4.1755, 12.6392, 35.8989, 7.0897, 18.0735,
    -20.1609, 19.4326, 6.9147, 47.0105, 43.7384, 47.9186, 42.4304, 34.0209, -25.9655, 19.7633,
    -22.5597, -0.5477, 27.7172, 52.3676, -41.2865, 12.1150, 13.5127, 9.0765, 59.9139, 23.5880,
    33.6844, 7.5006, 31.9522, 8.9824, -9.4438, -25.2637, -12.0464, 14.5995, 52.2297, 38.7223,
    25.2854, 44.4268, 55.7558, -1.9403, 17.3026, 13.9094, 13.4877, -13.8333, 43.9424, 0.3364,
    24.7136, 14.7167, 44.7866, -4.6796, 8.4657, 1.3521, 48.1486, 46.0569, -9.4456, 2.0469,
    -25.7479, 4.8517, 40.4168, 6.9271, 15.5007, 5.8520, -26.3054, 59.3293, 46.9481, 33.5138,
    25.0330, 38.5598, -6.7924, 13.7563, -8.5569, 6.1256, -21.1789, 10.6918, 11.8745, 39.9334,
    37.9601, -8.5211, 0.3476, 50.4501, 24.4539, 51.5074, 38.9072, -34.9011, 41.2995, -17.7334,
    41.9029, 10.4806, 21.0285, 15.5527, -15.4167, -17.8252
  ),
  longitude = c(
    69.2075, 19.8187, 3.0588, 1.5218, 13.2344, -61.8456, -58.3816, 44.4991, 149.1300, 16.3738,
    49.8671, -77.3963, 50.5577, 90.4125, -59.6162, 27.5615, 4.3517, -88.7590, 2.6289, 89.6390,
    -65.2627, 18.4131, 25.9231, -47.9292, 114.9400, 23.3219, -1.5247, 29.8739, 104.9282, 11.5021,
    -75.6972, -23.5087, 18.5582, 15.0445, -70.6693, 116.4074, -74.0721, 43.2551, 15.2662, 15.2663,
    -84.0907, -5.2893, 15.9819, -82.3666, 33.3823, 14.4378, 12.5683, 43.1456, -61.3870, -69.9312,
    -78.4678, 31.2357, -89.2182, 8.7832, 38.9251, 24.7536, 38.7469, 178.4419, 24.9384, 2.3522,
    9.4673, -16.5790, 44.8271, 13.4050, -0.1870, 23.7275, -61.7480, -90.5069, -13.5784, -15.5989,
    -58.1551, -72.3074, -87.2068, 19.0402, -21.9426, 77.2090, 106.8456, 51.3890, 44.3661, -6.2603,
    35.2137, 12.4964, -76.7931, 139.6503, 35.9284, 71.4704, 36.8219, 172.9790, 125.7625, 126.9780,
    47.9774, 74.5698, 102.6000, 24.1052, 35.4953, 27.4974, -10.7969, 13.1913, 9.5215, 25.2797,
    6.1296, 21.4314, 47.5079, 33.7738, 101.6869, 73.5093, -8.0029, 14.5146, 171.1845, -15.9582,
    57.5012, -99.1332, 158.1611, 28.8497, 7.4246, 106.9057, 19.2594, -6.8498, 32.5732, 96.1561,
    17.0832, 166.9315, 85.3240, 4.9041, 174.7633, -86.2362, 2.1254, 7.4951, 10.7522, 58.4059,
    73.0479, 134.4893, 35.2332, -79.5199, 147.1803, -57.5759, -77.0428, 120.9842, 21.0122, -9.1393,
    51.5310, 26.1025, 37.6173, -1.9536, -62.7177, -60.9789, -61.2248, -171.7667, 12.4578, 6.7273,
    46.7160, -17.4467, 20.4489, 55.4540, -13.2317, 103.8198, 17.1077, 14.5058, 159.9729, 45.3182,
    28.1879, 31.5825, -3.7038, 79.8612, 32.5599, -55.2038, 31.1367, 18.0686, 7.4474, 36.2765,
    121.5654, 68.7870, 35.7516, 100.5018, 125.5603, 1.2313, -175.1982, -61.5089, 10.1815, 32.8597,
    58.3260, 179.1942, 32.5825, 30.5234, 54.3773, -0.1278, -77.0369, -56.1645, 69.2401, 168.3273,
    12.4534, -66.8792, 105.8542, 44.2070, 28.2871, 31.0522
  ),
  stringsAsFactors = FALSE
)

# View the dataframe
head(capital_coords, 10)

# Check that Taiwan is included
capital_coords[capital_coords$iso2 == "TW", ]

# Optional: Sort by ISO2 code
capital_coords <- capital_coords[order(capital_coords$iso2), ]

# Print summary
cat("Total countries:", nrow(capital_coords), "\n")

library(plotly)
library(ggplot2)
library(maps)

# Create ggplot
world_map <- map_data("world")

p <- ggplot() +
  geom_polygon(data = world_map, 
               aes(x = long, y = lat, group = group),
               fill = "lightgray", 
               color = "white", 
               linewidth = 0.2) +
  geom_point(data = capital_coords,
             aes(x = longitude, y = latitude, 
                 text = paste0("Capital: ", capital, "\n",
                              "Country: ", iso2, "\n",
                              "Coordinates: ", round(latitude, 2), ", ", round(longitude, 2))),
             color = "red", 
             size = 2, 
             alpha = 0.7) +
  theme_minimal() +
  labs(title = "World Capital Cities - Interactive",
       x = "Longitude",
       y = "Latitude")

# Convert to interactive plotly
ggplotly(p, tooltip = "text") %>%
  layout(
    hoverlabel = list(bgcolor = "white", font = list(size = 12))
  )

now infer long and lats from either the succesful matches or if there is none…let’s just use the capital location…

  library(bigrquery)
  library(DBI)
  
  # Set your project ID
  project_id <- "patbis"
  
  # Option 1: Direct table download
  # Specify dataset and table name
  dataset <- "inglobe"
  table <- "person_address_extract"
  

  # Upload - this creates or replaces the table
bq_table_upload(
  x = bq_table(project_id, dataset, "lat_longs_mapbox"),
  values = lat_longs_mapbox,
  #fields = df,  # Optional: specify schema
  write_disposition = "WRITE_TRUNCATE"  # Replaces existing table
)



capital_coords=capital_coords %>% mutate(person_ctry_code=iso2) %>% select(person_ctry_code,latitude,longitude)
bq_table_upload(
  x = bq_table(project_id, dataset, "capital_coords"),
  values = capital_coords,
  #fields = df,  # Optional: specify schema
  write_disposition = "WRITE_TRUNCATE"  # Replaces existing table
)

  

anti=df %>% anti_join(lat_longs_mapbox) %>% select(person_id,person_ctry_code) 
anti=anti %>% filter(!is.na(person_ctry_code) & gsub("\\s+", "", person_ctry_code)!="")





lat_long=lat_longs_mapbox %>% select(longitude,latitude,person_ctry_code) %>% bind_rows(capital_coords) %>% group_by(person_ctry_code) %>% 
  mutate(rr=runif(1:n())) %>% arrange(person_ctry_code, rr) %>% mutate(nn=1:n()) %>% filter(nn<3) %>%  # Let's restrict to 5 places per country...
  select(longitude,latitude,person_ctry_code)


anti=anti %>% left_join(lat_long,by="person_ctry_code") %>%
  group_by(person_id) %>% mutate(rr=runif(1:n())) %>% arrange(person_id, rr) %>% mutate(nn=1:n()) %>% filter(nn==1)


anti=anti %>% select(-nn,-rr)

dfnew=anti %>% bind_rows(lat_longs_mapbox) %>% select(person_id, person_ctry_code,longitude,latitude)



write_parquet(dfnew,paste0(localbig,"\\data\\dfnew.parquet"))

library(leaflet)
library(dplyr)
library(arrow)
dfnew=read_parquet(paste0(localbig,"\\data\\dfnew.parquet"))


#df_sample <- dfnew %>% slice_sample(n = 10000)



library(data.table)
dfnew <- as.data.table(dfnew)
dt_sample <- dfnew[sample(.N, 100000)]

# With clustering (much better for 100k points!)
leaflet(dt_sample) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~longitude,
    lat = ~latitude,
    radius = 3,
    fillOpacity = 0.5,
    stroke = FALSE,
    clusterOptions = markerClusterOptions()
  )

bq_table_upload(
  x = bq_table(project_id, dataset, "dfnew"),
  values = dfnew,
  #fields = df,  # Optional: specify schema
  write_disposition = "WRITE_TRUNCATE"  # Replaces existing table
)